{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%reset -f\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from os import listdir\n",
    "from os.path import isfile, join\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "from scipy.integrate import odeint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def processData(fileName, varindex = 0, dataindex = 9, flgerr = 'stdev'):\n",
    "    working_dir = os.getcwd()\n",
    "    fileName = str(fileName)\n",
    "    callThis = join(working_dir, fileName)\n",
    "    \n",
    "    data = np.loadtxt(callThis, delimiter='\\t', unpack=True)\n",
    "    x = data[varindex]\n",
    "    y = data[dataindex]\n",
    "    x_sorted = np.sort(x)\n",
    "\n",
    "    idx_sorted = np.argsort(x)\n",
    "    y_sorted = y[idx_sorted]\n",
    "\n",
    "    var = []\n",
    "    avg = []\n",
    "    err = []\n",
    "    try:\n",
    "        for i, element in enumerate(x_sorted):\n",
    "            if i == 0: \n",
    "                i_begin = 0\n",
    "                prev_element = element\n",
    "\n",
    "            if element != prev_element or i == len(x_sorted)-1:\n",
    "            \n",
    "                i_end = i -1 \n",
    "                if i == len(x_sorted)-1:\n",
    "                    i_end = i\n",
    "\n",
    "                arr = y_sorted[i_begin:i_end+1]\n",
    "                if x_sorted[i_begin] != x_sorted[i_end] and i != (len(x_sorted)-1):\n",
    "                    raise Exception()\n",
    "                    \n",
    "                variable = x_sorted[i_begin]\n",
    "                var.append(variable)\n",
    "                avg.append(np.mean(arr))\n",
    "                if flgerr == 'stdev':\n",
    "                    e = np.std(arr)\n",
    "                else:\n",
    "                    ## this is standard error\n",
    "                    e = np.std(arr)/np.sqrt(len(arr))      \n",
    "                \n",
    "                err.append(e)    \n",
    "                \n",
    "                i_begin = i\n",
    "            prev_element = element\n",
    "            \n",
    "    except Exception:\n",
    "        print(\" --- something went wrong ----\")\n",
    "    \n",
    "    out = np.array([np.array(var), np.array(avg), np.array(err)])\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def twobody(t, initnum, gamma2):\n",
    "    return initnum/(gamma2 * t + 1)\n",
    "\n",
    "def onebody(t, amp, gamma):\n",
    "    return amp * np.exp(-gamma*t)\n",
    "\n",
    "def fullModel(n, t, gamma1, gamma2):\n",
    "    return -gamma1 * n - gamma2 * n**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "flgerr = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "two body loss rate = 0.694918348235 (1/sec)\n",
      "for 1288G, loss rate = 3.00913286594 (1/sec)\n",
      "for 980G, loss rate = 110.904526209 (1/sec)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAETCAYAAADKy1riAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4VGX2wPHvIQFCgCAdBJMgvUkLgiiIooiKbbGCYgFZUewoIro/dBe7oquIgAoKsYsNwQKKBUWKoIIuICWQBNDQIdTk/P54Z4aUSTIhycwkOZ/nuU+YO3fuPRlgzty3nFdUFWOMMSYQFUIdgDHGmNLDkoYxxpiAWdIwxhgTMEsaxhhjAmZJwxhjTMAsaRhjjAmYJQ1jjDEBC3rSEJERIrJERA6KyLR8jrtWRJaKyG4RSRaRJ0QkMoihGmOMySEUdxqpwH+AVws4Lhq4A6gDdAP6ACNLNjRjjDH5Cfo3d1WdCSAiCUDjfI6bmOVhiogkAmeUcHjGGGPyUZqae3oBK/09ISLDgGEAVatW7dKqVatgxmWMMaXe0qVL01S1bkHHlYqkISLXAwnAUH/Pq+pkYDJAQkKCLlmyJIjRGWNM6SciSYEcF/ZJQ0QuBh4DzlLVtFDHY4wx5VlYJw0R6QdMAc5X1d8KfYLevd3P+fOLMSpjjCm/gp40PMNmI4EIIEJEooAjqnokx3FnAonAJaq6KNhxGmOMyS0UdxoPAP+X5fHVwEMi8irwO9BGVTcCDwI1gNki4j32O1U9N5jBGlOWHT58mOTkZA4cOBDqUEyQREVF0bhxYypWrHhMrw/FkNuxwNg8nq6W5TgbXmtMCUtOTqZ69erEx8eT5cuZKaNUlW3btpGcnEyTJk2O6RxWRsSYcuzAgQPUrl3bEkY5ISLUrl27SHeWljSMKecsYZQvRf37tqRhjDEmYJY0jDEhdcMNN1CvXj3atWuXbf8999xDq1atOOmkk7jkkkvYuXMn4Drvr732Wtq3b0/r1q159NFHfa8ZP348bdu2pV27dlx11VW+ZpgjR45w//3307x5czp27EjHjh0ZN26c33hUlTPPPJPdu3cD7pv53Xff7Xv+qaeeYuzYsfn+TtOmTWPEiBG59p933nm+36Moxo4dy1NPPZXvMVdeeSVr1qwp8rVysqRhjAmp6667js8++yzX/rPPPpsVK1bw66+/0qJFC19yePfddzl48CC//fYbS5cuZdKkSWzYsIGUlBT++9//smTJElasWEFGRgZvvfUWAA888ACpqan89ttvLF++nO+++47Dhw/7jWf27Nl06NCBmJgYACpXrszMmTNJSyv63OLZs2dz3HHHFfk8gRg+fDhPPPFEsZ/XkoYxJmCJiYnEx8dToUIF4uPjSUxMLPI5e/XqRa1atXLt79u3L5GRboBn9+7dSU5OBtw3/3379nHkyBH2799PpUqVfB/w3n1HjhwhPT2d448/nvT0dKZMmcLzzz9PVFQUANWrV8/zbiExMZGLLrrI9zgyMpJhw4Yxfvz4XMd+8skndOvWjU6dOnHWWWexdevWfH/X+Ph4v8nnzTffpH379rRr145Ro0b59lerVo0xY8bQoUMHunfvnuv8a9eupXPnzr7Ha9asoUuXLgD07NmTuXPncuRItilwRWZJwxgTkMTERIYNG0ZSUhKqSlJSEsOGDSuWxFGQV199lXPPdVO0Lr30UqpWrUrDhg2JjY1l5MiR1KpVi0aNGjFy5EhiY2Np2LAhNWrUoG/fvvz555/ExsZSvXr1gK61YMEC3wev1y233EJiYiK7du3Ktv+0005j4cKFLFu2jCuvvPKYvtmnpqYyatQovvrqK5YvX87ixYv58MMPAdi3bx/du3fnl19+oVevXkyZMiXba5s2bUqNGjVYvnw5AFOnTuW6664DoEKFCjRr1oxffvml0DHlx5KGMSYgY8aMIT09Pdu+9PR0xowZU6LXHTduHJGRkQwaNAiARYsWERERQWpqKuvXr+fpp59m3bp17Nixg48++oj169eTmprKvn37mDFjRq7zTZ06lY4dO3LCCSewadOmXM9v3749V4KJiYlh8ODB/Pe//822Pzk5mXPOOYf27dvz5JNPsnKl30Lc+Vq8eDG9e/embt26vt/z22+/BaBSpUr0798fgC5durBhw4Zcrx86dChTp04lIyODt99+m4EDB/qeq1evHqmpqYWOKT+WNIwxAdm4cWOh9heH1157jVmzZpGYmOgbKvrGG2/Qr18/KlasSL169Tj11FNZsmQJc+fOpUmTJtStW5eKFSvyj3/8gx9++IFmzZqxceNG9uzZA8D111/P8uXLqVGjBhkZGbmuGRkZSWZmZq79d9xxB6+88gr79u3z7bv11lsZMWIEv/32G5MmTTqm+Q+qmudzFStW9P3eERERfpuaBgwYwJw5c5g1axZdunShdu3avucOHDhAlSpVCh1TfixpGGMCEhsbW6j9RfXZZ5/x+OOP8/HHHxMdHZ3tel999RWqyr59+1i4cCGtWrUiNjaWhQsXkp6ejqoyb948WrduTXR0NEOGDGHEiBG+D/WMjAwOHTrk97otW7Zk3bp1ufbXqlWLyy+/nFdeecW3b9euXTRq1AhwCe5YdOvWjW+++Ya0tDQyMjJ48803Of300wN+fVRUFOeccw7Dhw/n+uuvz/bc6tWradu27THFlRdLGsaYgIwbNy7bhzdAdHR0nkNXA3XVVVdxyimnsGrVKho3buz7UB4xYgR79uzh7LPPpmPHjtx0002A61/Yu3cv7dq1o2vXrlx//fWcdNJJdOvWjUsvvZTOnTvTvn17MjMzGTZsmC/2hg0b0q5dOzp16kTPnj259tprOf7443PFc/755zM/j8rYd999d7aO7LFjx3LZZZfRs2dP6tSpk+3YadOm0bhxY9/m7cjPqWHDhjz66KOcccYZdOjQgc6dO2friA/EoEGDEBH69u3r27d161aqVKlCw4YNC3WuAqlqmdq6dOmiPqef7jZjjF+///57oY6fMWOGxsXFqYhoXFyczpgxo4QiC53U1FQ966yzQh1GoTz55JP6wAMPZNv3zDPP6Msvv+z3eH9/78ASDeAzNqzX0zDGhJdBgwb5OqTLqoYNG3LjjTeye/du31DecHbJJZewdu1avvrqq2z7jzvuOK655ppiv54lDWOMyeHyyy8PdQgB++CDD/zuz9m/UVysT8MYY0zALGkYY4wJmCUNY4wxAbOkYYwxJmCWNIwxIeWvNHqoyqKbgpXppLF8uduMMeHLX2n0UJVFNwUr00nDGFO8EhMhPh4qVHA/i6PArb/S6KEqi24KZknDGBOQxEQYNgySkkDV/Rw2rHgSR36CWRbdFMyShjEmIGPGQI7K6KSnu/0lJdhl0U3BLGkYYwKSVwX0kqqMHoqy6KZgljSMMQHJqwJ6SVRGD1VZdFOwMps0EhMT2b17Nzt37Sy2tYyNKc/GjYMcldGJjnb7i8JfafRQlUU3AQikFG5xbsAIYAlwEJhWwLF3AluAXcCrQOWCzt+lSxedMWOGRkdH69d00K/poIBGR0eXyTLOxhRF4Uujq8bFqYq4n/ZfqnQqSml00XyWGiwJIvIPIBM4B6iiqtflcdw5wOvAmUAq8AGwUFXvy+/8CQkJmpaWRo+k4xnJASqQyXGkcT+x/BCX6neNXWPKqz/++IPWrVuHOgwTZP7+3kVkqaomFPTaoDdPqepMVf0Q2FbAodcCr6jqSlXdAfwbuC6Qa/RIOp4pLCeDCqyjDieQwhSW0yPJbkeNMaYowrlPoy3wS5bHvwD1RaR2Hsf7PMJGqrKfVBpxAd9QmSvpxH5S+ZkRI0bw+uuv+471t1C7McYY/8J5EaZquL4ML++fq5PjLkVEhgHDwI2uiCUVgLasoCJHiGAakMAuRjF9+nTWrVvH4MGDAbeI/MGDB2natClNmzalWbNm9OjRg969e5fk72aMMaVSOCeNvUDWtRa9f96T80BVnQxMBtensXFjBvGkEIUbYncXT/MY99ODbqxZcyrVqx/wvXbIkCGsXr2atWvXMmfOHLZs2cKNN95I7969ycjI8C0K700qTZs2pUePHrRq1aqkfm9jjAlb4Zw0VgIdgHc8jzsAW1W1oL4Q7ieWKWwnkwgAHmUMTVjPcF4gIUF4//0qdO3qOfb++7O9dt++fb7x3AcOHPCtv7tkyRLef/99jhw5wr///W8eeOABtm7dyhlnnJEtoTRt2pSEhATq1atXTG+DMcaEj6D3aYhIpIhEARFAhIhEiYi/5PU6MERE2ohITeABYFog1/ii9mpupCNHPDkxhXrMZyXVa5xHhQpw2mnw8sv+X1u1alVq167t+/OLL77I559/zp9//sn+/ftZu3YtQ4cOBVxSadmyJUlJSUyZMoXbbruN888/nzlz5gDw66+/0qdPH4YNG8bjjz/Oe++9x/Lly9m/f3/A75cxZd1zzz1Hu3btaNu2Lc8++ywAy5cvp3v37nTs2JGEhAQWLVoEuCkCt912G82aNeOkk07i559/9p1nzZo19O/fn6ZNm9KlSxfOOOMMvv3225D8TmVaIONyi3MDxgKaYxsLxOKapGKzHHsXsBXYDUwlwHkatWvXVkCXc6Iq6Cme69SuXVvT0lTPPlsVVIcMUU1PP+ahztlkZmbq5s2b9fvvv9ctW7aoquqPP/6o3bt317p162b7fefOnauqqvPnz9dBgwbpv/71L33ttdf0+++/182bN2tmZmbxBGVMAQo7T6O4/fbbb9q2bVvdt2+fHj58WPv06aOrV6/Ws88+W2fPnq2qqp9++qmefvrpvj/369dPMzMz9ccff9STTz5ZVVX379+vzZs3148++ijbuadOnRrsX6lUKMo8jaA3T6nqWE+S8KdajmOfAZ4p7DW2b98OQIbnRqp6lv21a8OcOfB//+dmsv78M7z3Hpx4YmGvkp2I0KBBAxo0aODb1717d3788UcAdu/ezdq1a1m7di2dOnUCYMuWLXz//fe8+eabZGZm+l63atUqWrRowaeffsrXX3+drekrLi7OVzLamKBLTHQVCjdudPVDxo0DTzHBY/HHH3/QvXt3X6mQ008/nQ8++AARYffu3QDs2rXLN3v7o48+YvDgwYgI3bt3Z+fOnWzevJnZs2dzyimncOGFF/rO3a5du2wLO5niUSY/fWJjY0lKSsqVNGI9RXIiIuA//4Hu3eGaa6BLF3j9dbjggpKLKSYmhk6dOvkSBsAVV1zBFVdcwaFDh9iwYYMvqcTHxwPuFv2FF17g4MGDvtdERkaye/duqlSpwsyZM1m/fr0voZx44olUrVq15H4JU755a6N7S916a6PDMSeOdu3aMWbMGLZt20aVKlWYPXs2CQkJPPvss5xzzjmMHDmSzMxMfvjhBwBSUlI44YQTfK9v3LgxKSkprFy5ks6dOxfp1zOBKZNJY9y4cQwbNoyMdNcRXg2Ijo7OtcRj//7uTmPAALjwQrj3XvfFKdhf5CtVqkSLFi1o0aJFtv1jxoxh9OjRpKam+hLK5s2bqVKlCgDvv/8+b7zxRrbXtG3blhUrVviezzqcuHbt2r5qocYUWn610Y8xabRu3ZpRo0Zx9tlnU61aNTp06EBkZCQTJ05k/PjxDBgwgHfeeYchQ4Ywd+5cb7N1Nv7+TV9yySWsWbOGFi1aMHPmzGOKzeQhkDas0rR16dJFVVVnzJihC2ivCvqvmjXzrTu1f7/qP//p+jl69lRNTs7z0LCzfft2Xbx4sb755pv6n//8Rx9++GHfc127ds3WlxITE6NXXXWV7/kPP/xQ582bpxs2bNAjR46EInwTYoXq0xBx/0lybiLFFs/o0aN1woQJGhMT4+vby8zM1OrVq6uq6rBhw/SNN97wHd+iRQtNTU3Vl19+WQcPHpztXIsXL/b1hZjsitKnEfIP+eLevElDVXV5TE/3K44bF9AbOWOGatWqqnXrqn7+eUAvCWvp6em6cuVK/fjjj3X8+PE6YsQIHZflvahfv74voVSqVElbtGihY8eO9T3/xRdf6O+//6779+8PRfgmCAqVNOLi/CeNuLgixbB161ZVVU1KStKWLVvq9u3btVWrVvr111+rqurcuXO1c+fOqqo6a9asbB3hXbt2VVX3b71p06bZOsK/+eYbSxp5KFUd4cGkUoFMhAp7cs0H9GvQIOjcGS6/HPr1g9Gj4aGHgt9cVVyqVKlCmzZtaNOmjd/nFy5c6Gv28m7VqrmxCAcOHOCcc85BVRERGjVqRNOmTRk6dChXX301GRkZLFu2jKZNm1KzZs1g/lomVMaNy96nAcVSG33AgAFs27aNihUrMmHCBGrWrMmUKVO4/fbbOXLkCFFRUUyePBmA8847j9mzZ9OsWTOio6OZOnUq4P6tz5o1i7vuuos77riD+vXrU716dR544IEixWZyK6Ufh4Hp2BFYEAF79wb8mtat4aef4Pbb4ZFH4Lvv4I03oHHjkoszVOLj44mPj6dPnz65nouMjOSHH37IllC8c1UAkpKS6OqZIVmrVi1fv8nw4cPp1asXBw4cYNu2bTRs2JAKFcK5xJkJmLffohhHTwF89913ufaddtppLF26NNd+EWHChAl+z9OqVStmz55dpFhMwcp00gDcUKkA7zS8oqNhyhTo3Rv++U/o0AGmTSvZ0VXhJjIyku7du9O9e3e/z9etW5cPPviAP//805dUFi1axGWXXQbATz/9RO/evYmKiuLEE0/0JZWbbrrJV+9LRKhUqVIwfy1TVIMGFTlJmNLNkkY+Bg2Crl3hiivc6Krbb4fHH4fKlYs5xlKoevXqXHzxxXk+37RpUyZMmOBLKOvWrWPevHlcfvnlALz77rtce+21nHDCCdmKRd5www3UqVPH1yxmjAkv5SNpFKJ5KqcWLWDhQjcc97nn4Jtv4K23oGXLYoyxDGrcuDE333xztn3ejjQ4Oj7fm1RmzpzJtm3bGDhwIACPPfYY48ePp1mzZtkmN15++eVUtqxdrCxBly/e/4PHKugr95W0hIQEXbJkiXvQuzf88ovrqPBMDiqKjz+GG26A/fvhv/91f7b/a8Vn165dVK9enQoVKjBnzhxmzpzpSyqbNm0iIiKC/fv3ExkZyT333MMXX3yRLaE0b97cb/9MKCUmJjJmzBg2btxIbGws48aNY1AYNe+sX7+e6tWr2xyeckJV2bZtG3v27KFJkybZngt05b6yf6cRGXnMzVM5XXihy0HXXANDh8Lnn8OkSWCDh4pHjRo1fH8+99xzOffcc32PDx48SEpKiq+EStOmTTnhhBP4448/mD17NgcPHuSEE05g48aNAAwfPpx169ZlSyqtWrUKakn7xMREhg0bRrpntFFSUhLDPDOowyVxNG7cmOTkZP7+++9Qh2KCJCoqisZFGNlT9u80/vc/iIqCYlwbPCMDnnwSHnwQGjSA6dPdpUxoZGZmkpqaSlpaGh07dgTgnnvu4euvv2bt2rXs3LkTgB49erBgwQLAraOSkZGRrfmrefPmxTp8OD4+nqSkpFz74+LibK16E3YCvdMo+0ljzRo4eBDS0or9WkuWwMCB8Oefrs/j4YfBBgOFn+3bt/Pnn3+SkZHBKaecAkD//v1ZtmwZqampvuOuuOIK3nrrLQCuv/56GjRokO1OpXHjxoUaPlyhQgW/7ccikq1ApTHhwJqnvIoweqogCQmwbBnceacbVfXFFzBjBuQxl86ESK1atTj55JOz7Zs1axYA6enprF+/nrVr1/rWUdm/fz8LFixgw4YNHD582PeakSNH8uSTT5Kens59992XLaE0adKEqKiobNfwFs7MyVs405jSqHwkjUOH3FYCtwFVq8LkyXD++a6fo0sXl0BGjACb0xb+oqOjadu2LW3btvXtq1KlCqtXryYjI4NNmzb5OuPbt28PQGpqKlOnTmVvllF5IsLkyZMZOnQomzdvZtq0aVxwwQW89NJLHDlyxHdcxYoVcxXONKY0KfvNU8nJsHata57yfJMsKVu2wJAhMHs2nHUWvPoqZKnibMoQVeXvv//ONlv+oosuolOnTnz55Zf07dvX7+sqVarE2LFj2bt3L40bN6ZRo0Y0atSIxo0bU7duXZs9b0LG+jTAJY3Nm2H1atcRHhdX4tdXdSOq7r4bKlaECRNcv4eNZixf9uzZQ8uWLdm8eXOu5+rUqcOOHTvIyMjItn/p0qV07tyZWbNmMX36dF8y8SaWrl272hwVU2KsT8PLW22whPo1chKBm25ydxrXXgtXXw0ffAATJ0LdukEJwYSB6tWrs2XLFr/Pbdu2jcOHD7N161ZSUlJISUkhOTmZpk2bApCWlsayZcuYNWuWb7guwNatW6lXrx5PP/00iYmJvoTi/Xn11VcTGRnJkSNHbHVHU2LK/p3G9u3w229ucp9n5EywZGTAU0/Bv/4FNWrASy/BP/4R1BBMCBV1yK2qsnPnTl9S6du3LxUqVOD111/n7bffJjk5mZSUFLZt20blypXZv38/IsINN9zAe++9ly2hNG3alAcffBCAjRs3EhUVRZ06daw5zPgEeqdR4L8YEYkSkS9EpHexRBZsEW71vqKUEinKpUeNgqVLXZXcAQPcnYdnCXNTxo0bN8639rWXvxUk8yIi1KxZk3bt2tGvXz/fB/zgwYP59NNP+eWXX0hLSyM9PZ1Vq1b5ZnSfe+65XH/99bRt25Y9e/Ywb948EhMTfecdOnQo9evXJyoqivj4eE477TTuvPNO3/PffPMNCxYsICkpiUOHDhX1bTBlTIH3sKp6QES6AhFBiKf4eZNGkJqn/GnXzpVbf+QRtzb53Lmu3+Oii0IWkgkC76zvki4jUqVKFeKy9NdddtllvmrDXllbFEaNGsUFF1yQrWksLcs8pptuuon//e9/vsf16tXj4osvZtKkSQBMnDiR6OjobHcy1atXL9bfyYSvgJqnROQ1YLeq3lryIRVNruap/fth0SJX2/zaa0MZGgDLl8N117lyJIMGuSKIJTyoy5hC+f3339m4caOv+SslJYUWLVowcuRIwPXX7M1x5z506FCmTJmCqvLPf/6TevXqZetzadKkCccdd1wofh0ToOLuCP8ceFJEGgKzga24ZUJ9VDU8Vz8JgzuNrDp2dDnskUfc+jVffgkvvuiarowJB/mt9giuQz41NdV3l5KSkkLr1q0BN1ly9uzZbN68Odus9/vvv59x48axc+dOLrjggmxDjRs1akS3bt2y3S2Z8BXonUZBNQ9UVcOi+SrXnUZmplt+r0kTWLcupLHl9MsvrlLuzz/DpZfCCy9A/fqhjsqYojty5IhvdFhycjLNmzenffv2JCcnc/XVV/v2HzhwAIAXX3yR4cOHs2LFCs4666xcc1gGDBjgW7zr8OHDvmWJTfEp7juNJgUfEqa8EyRyjIkPBx06uL6Op56CsWNh3jx45hnXimbzOkxpFhkZ6fvQz1rCpXHjxsyfPx9w/Sw7duwgJSWF+p5vS1FRUVx44YUkJyezYcMGvv/+e7Zv30779u1p2bIlc+fOpX///sTExGRLLKNHj6ZFixb89ddfpKam0qhRI+rUqWPl3ktA2R9yC7BggfsKn5wcsrgKsmqVK0Py/fdujsekSXDiiaGOypjQS09PJyIigsqVK/Pnn38yc+ZMX7OY9+fs2bM56aSTmDRpEjfddBPgZt97k0piYiKxsbH8+uuvrFq1yncH07BhQypWrBji37BkBbqmS7FP7hORysANQAJwAnCLqq4RkSuAX1X1j8B/jSCLiAjLO42sWrZ0qwK+9JIbptuunauae8cdR+cnGlMeZR223KxZM+699948jz3vvPN4//33syWU5ORkX3PWu+++y3/+8x/f8SJC/fr1WbVqFTExMXzyySesWLEi12z80tocVhJrugTap9EC+BKoASwFegNdVfVnEXkBiFHVwQFdUKQW8ArQF0gDRqvqG36Oqww8B1wCVAQWADepakp+5892p+FVtSpER0MpWWhm0ya45Rb45BPo3BmmTHE/jTFFs3v3bpKSkrIllC1btjBx4kREhOHDh/PSSy9le010dDR79+5FRHjqqaf4/fffsyWUuLg4XzHLcFOYCabFWntKRD4DqgIXAHuBQ0CCJ2lcBjyuqgE1pojIm7hJhUOAjsCnQA9VXZnjuHuBQbjksguYAlRV1XznVPtNGjEx7ut6KZpVpwrvvQe33upy3R13wEMPQSn9wmNMqZGenp5tDsu+ffv45z//CcCIESP48MMPs40Oa9GiBatWrQLcxMsNGzZku1Np06aNr4BlZmZmUGfhF2ZNl+JOGvuAy1R1tohEAIc5mjR6AZ+rapUAzlMV2AG0U9XVnn3TgRRVvS/HsROBPap6r+fx+cAzqtoyv2v4TRo1a7pRVLt2Ffi7hpudO2H0aNdsdcIJboTVhReGOipjyrcjR46wZcsWUlJSOHjwIL169QLg3nvvZdGiRb47mQMHDtCnTx/mzp0LQKtWrdi9e3e2kWGnnXYaV155JeDWbK9Xrx5Vq1YtljhL4k4DVS1wA7YBAzx/jgAygc6ex1cBmwM8Tydgf459I4FP/BybgGuSOh6IBt4Ans3jvMOAJcCS2NhYzaV2bdWqVXPvL0UWLFBt104VVC+6SDUpKdQRGWPyk5mZqWlpabpp0ybfvkceeUSHDBmi/fr103bt2mnNmjX12muv9R0fFRWlgB533HHatm1b7du3r7766qu+5z/99FNdvny5/v3335qZmVlgDDNmzNDo6GjFzatTQKOjo3XGjBm5jgWWaACf44F2sX4J3C8ic3HNUwDq6Xe4FTfhLxDVcE1NWe0C/NUgWA1sBFKADOA3YIS/k6rqZGAyuDuNXAeUgo7wgvTo4eZzjB/vmqlat4b/+z+3amAZH/xhTKkkIr7VIL1Gjx6d6zhvifzMzEymTJmSqxN/27ZtgKuOfP755/teV7lyZRo1asT999/PkCFD2L17N1OnTs3WNHbFFVcAxVvKJtCkcQ/uW/+fuASiwL+AtkAlINDarXuBmBz7YgB/07UnAlFAbWAfcC8wB+gW4LWOiows9UkDXHK491644gq47TY3yur1192aHaefHurojDHHIsJTtSIiIoKrr746z+NiYmJsT9wJAAAgAElEQVRYuHAhycnJ2ZJKvXr1AFi3bh133HFHtteIiG+0lKqybds23nnnHdLS0qhXrx716tWjU6dO1KpVK+B4A56nISI1gbuAPkAdYDswD9fPsC3Ac3j7NNqq6hrPvteBVM3dp7ECGKOqH3keH+d5bV1VTSMPfvs0YmPdHI3Mgia2ly4ff+ySR1KSq5775JPQoEGoozLGhIKqsn379mwJ5bPPPmPOnDm+mff+fPbZZ5xzzjnhu3KfiLyFu1MZihs9NRv/o6em4u5CbgDScXc7t6hqo/zO7zdpNGniVu47dKjMteWkp7s6Vk88AVWquLkdt9xiczuMMXl3hDdu3Jgvv/ySv/76i/bt21OzZs3iW08jKxE5TkROE5HLRORUz7f/wroZqAL8BbwJDFfVlSLSU0Syls4cCRwA1gB/A+fh5mwUXgjX1Chp0dGu3PqKFW6NqTvugE6d3ERBY0z5tnHjRr/7U1JSaNWqFb169aJmzZqFOmdASUNEIkXkcSAZ+BZ4G/gOSBaRJ0Qk4K/vqrpdVS9W1aqqGqueiX2q+p2qVsty3DZVHaSq9VT1OFU9TVUXFeq38wqzSrcloUULmDPHLS27Z4+roHLVVW6ioDGmfIqNjS3U/kAEeqfxDHA78AjQBten0QZ4FLgNePqYIyhpiYlHq9t27+4el1EicPHF8PvvbmTVhx9Cq1auBHs+TZrGmDKqqKtH+hXIuFxcB/RdeTx3N7AjkPMEY+vSpUvWQcqq0dFucoN3i452+8uB9etVBwxwv3Z8vOp776kGMLTbGFOGzJgxQ+Pi4lRENC4uzu8cDdXA52kEOiN8GzBQVT/389w5wJuqGviYrRKUrSM8Pt4NLcopLs51jJcT8+a5vo4VK1yz1bPPurLsxhjjVdwd4dNxo538uRGYEWhgQZVHJ1Ce+8uoPn1g2TK3QuBvv7mO8mHDYOvWUEdmjClt8kwaInKzdwM2AKeIyEoReVRE7vT8/B032W5tkOItnLw6e4rQCVRaRUbC8OGwZo2765g6FZo3h8cey7+/o3fvo8uSGGNMns1TASzxmpVqOC73mpjovlJ7askDbozq5MlQhGn0ZcHq1TBypCu/HhcHjz4KV16Ze8VAb8LwLLZmjCmjitw8paoVCrGFRcLIZdAglyAqVXKPIyIsYXi0aOFmlM+b54oADxzoBpd9912oIzPGhLPgFXYPlUGD3KchQKNGljByOPNMWLLENVelpECvXnDJJW75WWOMyalQxSZEpCXQCFdIMBtVDbTSbfB521zKQNHCkhARAdddB5df7kZWPfYYtG3rmqZ+/NFVX4mPd/M9LOcaU74FlDREpD2u5EdrQPwcorh1NsKXiCWNAkRHw/33w9ChrgDil18efS4pyXUPgSUOY8qzQJunXsWt1tcfaAk0ybEFtNRrSFnSCFi9eq6jPKf0dJdUjDHlV6DNU61xK/flmtxXaljSKJT8prhMn+46ziPC+97SGFMCAr3TWASU/skNR46EOoJSI6+pLJUqweDB0LEjfPSRq81ijCk/Ak0aw4BhIjJIRI4XkeicW0kGWSzsTqNQxo1zfRxZRUfDyy/D22+7zvGLL3bL0M6bF5oYjTHBF2jSSMPNCn8d2IRbnjXnFt4saRSKd4pL5crucVyce3zNNW6U1cqVMGWKWxDxrLPc0N0ffwxtzMaYkhdon8YM4BTgKdw64YdKLKKSYkmj0AYNcokBcs8Ij4w8Ospq0iS3emCPHnDuuW71wIQC55UaY0qjQJPGGcCN6lkwqdSypFHsoqLg9ttdAnn+ebdOedeucNFF8NBDVk3XmLIm0OapDbh1uksvEcjMtM7wElK1Ktx3H6xfD//+t7sz6dgRBgyAX38NdXTGmOISaNK4BxgjIvElF0oJ884K37cvtHGUcTEx8MADbrmS//s/mDvX3W0MGAC//BLq6IwxRRVo0ngIN+R2tYisFpFFObcSjLF4eJNGGV4nvCTMn39sFW6POw7GjnXJ41//csmjY0dX1+rnn4s3RmNM8ASaNFYAs4FEYAGw0s8WvubPd8WTwJJGkNWs6fo2kpJcEpk/H7p0gf794aefQh2dMaawAlrutTTJtp6GV2Kiq8h35Ag0aABPPWUFlEJk1y7XYT5+PGzfDmef7ZqzevUKdWTGlG/Fvdxr6eVdiMnbAb5li5vSnJgY2rjKqRo1jvZ5PP646+c4/XTo2RM++8xmmBsT7gK60xCRdwo6RlUvL5aIiijXnUZ8vGsbySkuzn1ymZBKT4dXXoEnnnATBTt3htGjXd+H1bYyJniK+06jrp+tJXAhcCpQ5xjjLHn5Vd4zIRcdDbfeCmvXuhIle/bAZZdBmzbw6qtw8GCoIzTGZBVQ0lDVM/xsHYDmwGZgfIlGWRR5Vd7La78JiUqVYMgQ+OMPeOcdN+9jyBA48UTXBbV7d6gjNMZAEfs0VHUT8CjwRKCvEZFaIvKBiOwTkSQRGZjPsZ1F5FsR2SsiW0Xk9kIH6a/yXoUKbr8JOxER7k5j6VL4/HNo1Qruucfl+Pvvd11SxpjQKY6O8AygcSGOn4CrXVUfGARMFJG2OQ8SkTrAZ8AkoDbQDPii0NHlrLxXoQK0aGGjp8KcCPTt6yroLlrkiiI+9pjriho61N2RGGOCL9CO8DZ+dlfCLc70b2Cjqp4ZwHmqAjuAdqq62rNvOpCiqvflOPYR4ARVvabAALPwO+QW3ILXy5a5pNGhw7HNWDMh9eef8MwzMHUqHDjg5nqMHOmG64q/RYiNMQEr7o7wFcBvObaluMl+24ChAZ6nBZDhTRgevwC57jSA7sB2EflBRP4SkU9ExG9HhIgME5ElIrLk77//zvvqFSvC4cMBhmrCTbNm8OKLbgzD2LGwcKH7LtC1K7zxhv3VGhMMgSaNM4Azc2w9cHcC3VR1XYDnqQbsyrFvF1Ddz7GNgWuB23ElTNYDb/o7qapOVtUEVU2oW7du3le3pFEm1K3r6lpt3AgvvQR797rWxhNPdEN3d+wIdYTGlF2Bjp76xs/2k6qmFPJ6e4GYHPti8L+I037gA1VdrKoHcPWveohIjUJe8yhv0rAZZGVClSrwz3/C77/DrFnQvDmMGgWNG8Mtt8Dq1QWfwxhTOIXqCBeRyiJyooi0ybkFeIrVQKSINM+yrwP+a1f9CmT9dPf++dhbrytWdAkjM/OYT2HCT4UKcP758NVXsHy5W1nw5ZehZUu3//PP7XuCMcUloKQhIo1EZBZuTY01ZO/b8PZ3FEhV9wEzgYdFpKqInApcBEz3c/hU4BIR6SgiFYEHge9VdWcg1/KrYkX305qoyqwOHVxHubffY+lS6NfPTRZ88UXXlGWMOXaB3mlMARKAu4B+ZO/b8PZ3BOpmoArwF66PYriqrhSRniLi+y+tql8B9wOfeo5tBuQ5pyMg3k+Mn35y5UWs/lSZVb++6/dISoLp06FaNddk1agR3HEHrFkT6giNKZ0CHXK7C7fca4E1qEItzyG3rVvDqlXZ2ymio90cDpuzUeapuu8Kzz8P777rbjb79XOJ5Nxzrc6VMcU95PYvXMd06bV+fe6G7fR0GDMmNPGYoBKB7t3dzeXGjfDww24Z2gsucEN5H38c8hutbYxxAk0a/wJGiUjOkU+lR16V76xwYbnToAE8+KArcvzOO66l8r773Kira66BH36wjnNj8hJo89S7QDfcfIrFQM7OaFXVK4o/vMLLs3kqKsp/4rAS6QY3bHfiRHjtNVdp96ST4KabXMtlTOn9qmRMwIq7eaoOsBZYDlQkd5n0escYZ/A0aeLGZmYVHW2FCw3gRlc9/zykprpurogIuPlmOP54uPFG8Pc9xJjyqHws9wqu3sTWrUc7w+PiXMKwTnDjhyosXgyTJsFbb7nur06d3CKQAwfa3Ycpe2y5V3/q13f/22vUcE1SljBMHkTg5JPdqoKpqTBhgpsTOnw4NGwIN9xgfR+mfCpfSQOs/pQptBo1XFPVsmWuTPvAga4D/dRToW1bV3n3r79CHaUxwWFJw5gAibiKulOmwObNrlRJjRpw991u0uCAATB7Nhw5EupIjSk55TdpWLuCKYLq1d1ytD/+CCtWwG23wbffulpXcXEwerTrPjOmrCk/SWP+fLd560/tylmh3Zhj07YtPP00pKTA++9D587w5JNuqdoePdxorJ3HXjHNmLBSfpKGlzdppKWFNg5T5lSqBP/4B3zyCWza5Nb22LXLlW9v2BCuugrmzLHmK1O6Reb1hIi8A4xW1bWeP+cnbCb3FSjS8yunpbn6EcaUgIYN4Z573HK0S5fCtGnw5ptu+G6DBm7g3jXXuKq8xpQm+d1p1MVN5AM3eS/nhL7SNbnPy3unsW1baOMw5YIIJCTACy+4zvOZM10NrOeeg44d3czzJ590TVvGlAZ53mmo6hlZ/tw7KNEEgzVPmRCpVAkuucRtaWnw9tswYwbce69bcfCMM9wdyIABblSWMeGoyH0aItJLRL4qjmCCwpKGCQN16riy7D/+6Jal9a55PmSIm4N66aXuruTAgVBHakx2xdERXhc4vRjOExwREa7NwJKGCRPNm7uksXo1LFzoOs6/+87dcdSvD9dfD198YR3oJjyUv9FTIu5uw5KGCTMi0K2b6+9ISXFrm19yibvjOOccN4FwxAj4/ntb5t6ETvlLGvPnQ8uWljRMWIuMhL593airrVvd/I9evVwtrJ493QTCu+5yqxHaPFUTTOUvaQBkZMBHH7mvdrZWuAlzUVFu/se777oaV4mJruLuhAluJFaTJm5476JFZTOB9O7tNhMe8puncXOA5+hYTLEER2Ji9rXCk5JcvWuwqrcm7FWv7gomDhzoZpl//LEbhfXcc/DUU+470KWXwmWXuTpZIqGO2JQ1ea6nISKFaTVVVY0onpCKJs/1NLzi412iyMlW8DMlIDHRLUO/cSPExpbcEi47drib53ffPdppHhvrOtMHDIBTTsm9Bllp4b3LmD8/lFGUfYGup1F+FmHyqlDB/z28iPUummKVmOhuYtPTj+6Ljna1qErypnbHDlfK5L33XGf6oUNuFvoll7gE0qvX0ZHnpYEljeCwRZjyEhtbuP3GHKMxY7InDHCPx4wp2evWrAmDB7umq7//hjfecGt/vPYanHWWSyDXX++e37+/ZGMxZU+hkoaIRIrIiSLSJudWUgEWu3Hj3NTcrGytcFMCNm4s3P6SEBPjCiW+955LIB984Mq3f/ABXHSRm2Q4YABMnw7btwcvLlN6BZQ0RKSiiEwEdgNrgN/8bKXDoEFw551HH8fFlXx7gSmXwu2mNjoaLr4YXn/djcL64gu47jo3oXDwYKhXD848E/77X+veM3kL9E7jX0B/YAggwAjgemAesAG4oCSCKzGXX+5+fvihrRVuSsy4ce6DOqtwuamtVAnOPtsN2920yc33GDXKzQm5/XY3jLdDB3jwQVi82Lr7zFGBJo3LgbGAt0T6IlV9XVX7At8DFwV6QRGpJSIfiMg+EUkSkYEFHF9JRP4nIsmBXqNAdeq4nzbBz5SgQYPcTWxcnBtnEa43tRUqwMknu2S2cqUrZ/L0065v5JFH3HONGsGNN7p+kH37ghdbYqK7E/rmG5tSFS4CTRonAKtVNQM4ANTM8lwiMKAQ15wAHALqA4OAiSLSNp/j7wH+KsT5C+ZNGn//XaynNSanQYPczWxmZum5qW3e3M02nz/fNWNNn+5GXL39tusHqV0bzjsPXnzR/+j14uIdfXbwoHvsnVJliSO0Ak0am4HjPH9eD/TK8lzTQC8mIlVxCeZBVd2rqt8DHwPX5HF8E+Bq4NFArxGQ6Gg3Syo1tVhPa0xZU7s2XH21SxhpaTB3Lgwf7u5GbrnFfftv1841bX37LRw+XHzXDtXoM5O/gOZpiMgrwDZVvVdE7gCewjVVHQSuAN5U1SEBnKcT8IOqVsmybyRwuqrm6hcRkVnAK8AOYIaqNs7jvMOAYQCxsbFdkgL5+tO+vVsVp107GwBuTCGpusTx6adu+/ZbN6GwRg1XM+vcc6FfP7eC4bGyKVXBFeg8jTzLiOQwBqgDoKrPiogAlwJVgOeBhwM8TzVgV459u4DqOQ8UkUuASFX9QER653dSVZ0MTAY3uS+gSOLjYf36gA41xmQn4up+tmzpmrJ273Z3IbNnu+3dd91xHTu6BHLuua5OVmEmFcbG+m/+silVoRVQ85SqblHVFVkej1fVU1W1s6qOUtVAu8b2AjE59sUAe7Lu8DRjPQHcGuB5Cy8u7mhjqTGmSGJiXFHFl192Zd2XLXOd6NWrwxNPuD4R75yQyZMD6wsJ59Fn5VmgdxrFZTUQKSLNVXWNZ18HYGWO45oD8cB37qaGSkANEdkCdFfVDUWOJC7O3U/byjbGFCsRd4fRsSOMHg27dsFXX8Fnn7lt5kx3XIsWbp2Qvn1dqZBq1bKfxztoYMgQ9/0uLq7kaneZwOVXsLAwS7iqqvYJ6IIibwEKDMVVyJ0N9FDVlVmOicTTHObRA3gB6Az87RnF5VeBtae83nkHrrgCunSBQI43xhSZqisy/fnnbps/35UyqVjRFVU8+2y3JSS4RTbBak8FS3H0aWwL4DoNcR/ohal6eDPwKm4Y7TZguKquFJGewBxVraaqR4At3heIyHYgU1W3+D3jsYiPdz9tEWZjgkYEWrVy2+23uzuIBQvc7PQvv3STCR98EI47zs1OP+ssN2KqSpWCz22C45iq3IpILDAKuAHXHzFeVYt3WOwxCvhOY+tWV7mtWTNYs6bg440xJe7vv2HePNep/uWXR+t0Va7samj16eOSyfHHhzbOsqhESqOLSDNgNG7uxF/A08AkVQ2bWpkBJw1Vt6bm8ce7OgrGmLCiCmvXuomEO3a4YbbeooqtWrnkceaZrvmqdu2QhlomFOuQW8+M7THAZcAm4HbgVVU9VKQoQ0nEfX2x5iljwpKIawg4/ni3ffUVLF/ufn71lSv1/uKL7tgOHeCMM9zWs6crgWJKRr5JQ0S64JLFRbiRT0Nxk+zy7IguVaKibNitMaVEhQrQubPbRo50s88XL4avv3ZJ5KWX4Nlnj47eOv10dxfSsyfUqhXq6MuOPOdpiMgcYBHQBLhSVVur6mtlJmEkJrqxgHv2WCU0Y0qhihWhRw9XVmTePLdm+jffwNixbmb6xImuFHydOu5O5Lbb4P33reRcUQWyRvh2oMBJ+6parxjjOmYB9WmEah1OY0yhHeuQ24MHYdEi97pvvoEffji6UmHr1m7CYc+ebrNZ5sXQES4i/1eYC6rqQ4U5vqQElDTi4/1PSY2Ls9VnjCmjDh1yU7K++87Vyvr+e1f+BFzS6NkTTjvNbW3auOaw8qRERk+VBgElDauEZky5l5EBv/3mkoh32+KZCVazpmv6Ou00t756166uC7QsK+6ChWWLVUIzptyLiDha7uTWW933yHXr3B3Id9+5SYeffuqOrVjRFY/o0ePoVpQKvqVZ+bzTsD4NY0wA0tJcX8iCBe7n4sVHB1zGx7vSJ6ec4pLISScVropvuLHmqYIkJh6thGYJwxgTgIMHXQXfH390SeSHH46u5ValimvG6t796Faa7kYsaQSid283W6hdO3dPaowxhaDqCkr8+KPbFi6En38+uoLhCSdAt25Hty5dcpd7DxfWpxGoqCgbMWWMOSYiris0NtYVzQZXZGLZMvjpJ5dEfvoJ3nvPPRcR4b6jdusGJ5/stjZtjlb0LQ0saURFuapohw5BpUqhjsYYU8pFRR3t6/DautXNGVm0yCWRd95xLeIAVau6We5dux7dTjzRJaTikJjoJkBu3OiSW1HXJLGkUbmyu8dMTnZ/U8YYU8zq14cLLnAbuJH9f/7pksjixe7nhAlHO9lr1nRriiQkuCathAT3gV/YRJJzzE9SknsMx544rE9jxw749VdXh+DMM0s0NmOMycvhw7BihZuAuHix21asOLq4aO3aRxNIly7u7iQuLv9EUph5zNanESjvjJ1AFi02xpgSUrEidOrkthtvdPsOHHATEBcvdh3sS5e6Nde9iaRWLXe8t5Bjp07QvPnR2eze9Uhyymt/ICxpVK7sfj70EFx/fWhjMcaYLKKijvZzeHkTiTeJLF0Kzz3numXBrbXeoYNLILVqwTY/a7AWZR6zJY0KFVwHuK2rYYwpBfwlkkOH4I8/XCJZtsxt06bB3r25Xx8d7TrDj5UlDXB/C5Y0jDGlVKVK7u6iQ4ejDSaZmW7lw+efdwtW7d7t+jKKOnqqfHeEex1/vCuyf+hQ8Y1zM8aYUiTQjvByVvw3D9WquZ4l6ww3xph8WdIAlzTANQQaY4zJkyUNcFMywZKGMcYUwJIGuMIv0dGWNIwxpgCWNLyqVbOkYYwxBbCk4VWtGqSkuFFUxhhj/Ap60hCRWiLygYjsE5EkERmYx3H3iMgKEdkjIutF5J4SDcw6w40xpkChuNOYABwC6gODgIki0tbPcQIMBmoC/YARInJliUVlScMYYwoU1KQhIlWBAcCDqrpXVb8HPgauyXmsqj6hqj+r6hFVXQV8BJxa7EElJrqVUn74wXWIf/hhsV/CGGPKimDfabQAMlR1dZZ9vwD+7jR8RESAnsDKYo3GW2zeW8Q+I8OtkJKYWKyXMcaYsiLYSaMasCvHvl1A9QJeNxYX61R/T4rIMBFZIiJL/i5MR/aYMUdXJ/FShdGjAz+HMcaUI8FOGnuBmBz7YoA9eb1AREbg+jbOV9WD/o5R1cmqmqCqCXXr1g08mryKym/aFPg5jDGmHAl20lgNRIpI8yz7OpBHs5OI3ADcB/RR1eRijyavovK1ahX7pYwxpiwIatJQ1X3ATOBhEakqIqcCFwHTcx4rIoOAR4CzVXVdiQQ0bpybCZ5T+/YlcjljjCntQjHk9magCvAX8CYwXFVXikhPEcm6ZMh/gNrAYhHZ69leKtZIBg2CyZOPrt4XFwft2sGePFvLjDGmXLP1NAB693Y/58+HUaPg2Wdh506oUqW4wzPGmLBk62kcqzPPdIsxff11qCMxxpiwY0kjp9693ezwTz4JdSTGGBN2LGnkVLky9O0Ls2a5ORvGGGN8LGn4078/JCfD8uWhjsQYY8KKJQ1/zj8fRKyJyhhjcrCk4U+9etCtm2uiMsYY42NJIy/9+8PixbB5c6gjMcaYsGFJIy8XXOB+fvppaOMwxpgwYkkjL+3bu9pU1q9hjDE+ljTyIuLuNr78EvbvD3U0xhgTFixpgCsfMn9+7v0XXOAShnWIG2MMYEkjf336QNOm8OSTNtHPGGOwpJG/yEi49143ispqURljjCWNAg0eDA0awGOPhToSY4wJOUsaBYmKgjvvdB3iS5eGOhpjjAkpSxqBuOkmqFHD7jaMMeWeJY1AxMTALbfA++/D6tWhjsYYY0LGkkagbr/dlU2/6y4bSWWMKbcsaQSqXj144glXVuS550IdjTHGhIQljcIYMQIuusgNw7VOcWNMOWRJozBE4NVX3RDcK66A3btDHZExxgSVJY3CqlUL3ngD1q93czgOHgx1RMYYEzSWNI7FaafB+PHw0UduPfHt20MdkTHGBIUljWN1223ujmPhQjj1VHfnYYwxZZwljaK46iqYOxe2boWuXWHiRDhyJNRRGWNMibGkUVQ9e7q7jbZt4eab3eJNn3xiczmMMWWSJY3i0KKFW4/jww8hMxMuvNCVVH/gAfjjj1BHZ4wxxSboSUNEaonIByKyT0SSRGRgHseJiDwuIts82xMiIsGON2Aibg7HihXw2mvQrBk8+ii0aeMSyDXXuOarxYth585QR2uMMcckMgTXnAAcAuoDHYFPReQXVV2Z47hhwMVAB0CBL4F1wEtBjLXwKlZ0Q3EHD4YtW+Ddd91dyNy5MGPG0ePq1HHJpGFDN++jfn2oWdMVRqxRA6pVg+hoqFLFVdqtXBkqVXLnr1jRrfURGQkVKkBExNGfYZxXjTGln2gQ295FpCqwA2inqqs9+6YDKap6X45jfwCmqepkz+MhwI2q2j2/ayQkJOiSJUtKJP4iUYUNG+DXX2HNGlf4cP16l1i2bIG0tOK7loj/LetzWY/N+dq8nsvvesaY8Na/P7z5Zp5Pi8hSVU0o6DTBThqdgB9UtUqWfSOB01X1ghzH7gL6qupPnscJwNeqWt3PeYfh7kwAWgKrsjxdByjGT+QSEe4xWnxFF+4xWnxFF+4xFhRfnKrWLegkwW6eqgbsyrFvF5ArEfg5dhdQTUREc2Q6z93IZH8XFJElgWTPUAr3GC2+ogv3GC2+ogv3GIsrvmB3hO8FYnLsiwH2BHBsDLA3Z8IwxhgTPMFOGquBSBFpnmVfByBnJziefR0COM4YY0yQBDVpqOo+YCbwsIhUFZFTgYuA6X4Ofx24S0QaicjxwN3AtGO4rN9mqzAT7jFafEUX7jFafEUX7jEWS3xB7QgHN08DeBU4G9gG3Keqb4hIT2COqlbzHCfA48BQz0tfBkZZ85QxxoRO0JOGMcaY0svKiBhjjAmYJQ1jjDEBK7NJI9AaV0GOab6IHBCRvZ5tVZbnBnri3CciH3r6fko6nhEiskREDorItBzP9RGR/4lIuoh8LSJxWZ6rLCKvishuEdkiIncFMz4RiRcRzfI+7hWRB0MQX2URecXz97ZHRJaJyLlZng/pe5hffOHyHnquNUNENnuutVpEhmZ5Lhz+HfqNL5zeQ8/1mns+X2Zk2Zfn58oxf0aqapncgDeBt3GTBE/DTQ5sG+KY5gND/exvi5ur0ssT7xvAW0GI5x+4+l4TcSVbvPvreN6vy4Ao4ElgYZbnHwW+A2oCrYEtQL8gxhePq0cWmcfrghVfVWCsJ54KQH/P32N8OLyHBcQXFu9hln//lT1/buW5VpdweA8LiC9s3kPP9b7wXG9Glrjz/FzhGD8jSyT4UG+e/yyHgBZZ9k0HHgtxXPPxnzQeAd7I8ripJ/7qQYrrPwrr+nkAAAmlSURBVGT/UB6GK/eS9f3cD7TyPE7BlXjxPv9vSjDJ+YmvoP+sQY0vx7V/BQaE23voJ76wfA9xZYA2A5eH43uYI76weQ+BK4F3cF8SvEkjz8+VonxGltXmqRZAhnqKInr8gsu8ofaoiKSJyAIR6e3Z1xYXHwCquhbPX2gI4vMXzz5gLdBWRGoCx2d9ntC9t0kikiwiU0WkDkAo4xOR+ri/s5WE4XuYIz6vsHgPReRFEUkH/of7UJ5NGL2HecTnFdL3UERigIdxc9myyu9z5Zg/I8tq0ihMjatgGgWcCDTCTbT5RESaEn7x5hdPtSyPcz4XLGlAVyAO10xQHUj0PBeS+ESkoieG11T1f4TZe+gnvrB6D1X1Zs/5e+ImAB8kjN7DPOILl/fw38Arqropx/6C3r9j+swpq0mjMDWugkZVf1LVPap6UFVfAxYA5xF+8eYXz94sj3M+FxSquldVl6jqEVXdCowA+nq+cQU9PhGpgLu1P+SJBcLoPfQXX7i9h56YMlT1e6AxMJwweg/9xRcO76GIdATOAsb7ebqg9++YPnPKatIoTI2rUFJAyFFnS0ROBCrjfo9QyBlPVVx76EpV3YG7PQ+numDeGaoS7PhERIBXcIuKDVDVw56nwuI9zCe+nEL2HvoRiee9Igzew3ziyykU72FvXN/KRhHZAowEBojIz+T/uXLsn5El2WkUyg14Czc6oCpwKiEePQUcB5yDGwUSCQwC9uE61toCu3G3vlWBGQSn0zHSE8+juG+i3tjqet6vAZ59j5N91MpjwDe4USGtcP85SmLUSl7xdfO8bxWA2rgRIF8HOz7PtV4CFgLVcuwPl/cwr/jC4j0E6uE6casBEZ7/I/twNelC/h4WEF/I30MgGmiQZXsKeM/z3uX7ucIxfkYW+3+icNmAWsCHnr/gjcDAEMdTF1iMu/3b6fmPfHaW5wd64twHfATUCkJMY3HfjrJuYz3PnYXr9NuPG/UVn+V1lXH1w3YDW4G7ghkfcBWw3vNebcYVt2wQgvjiPDEdwN3ue7dB4fAe5hdfGL2HdXEfrDs91/oNt0Kn9/lQv4d5xhcu76Gf/zMzsjzO83OFY/yMtNpTxhhjAlZW+zSMMcaUAEsaxhhjAmZJwxhjTMAsaRhjjAmYJQ1jjDEBs6RhjDEmYJY0TFCJyFgRScvjuWkisqSQ5/OuadC/gONGiEi+48tFpK+I3FEccZUEEYn2rMtweqhjCYSIVBGRv0SkZ6hjMcXHkoYxR/UFciUNXEG464Ibil+3AutV9ZtQBxIIVd0PPI97/0wZYUnDmAKo6lpVXRHKGDxFB2/BzTAuTaYBvUSkfagDMcXDkoYJayLSUUTmeZb73CEiiZ51IfJ7TWUReUFEdorIdhEZD1Qs4DVjcesRxHmau1Q8S8zmbJ4Skes8z3cWt4Rvuogs9zyu6llXYZeIrBORq/xc6yJxy9ge8DQ3PeEpXZ6fM3El9WfmOJeKyJ0i8rSIbPOs1TLS89y1nhh2ilt2NCrL644TkZdFJNUTx0YRmZLj3O1E5FNxS8XuEZF3RaRBjmNqi8gkccuhHhCRVVmb+NSV614MDC7g9zOlRGSoAzDlk4j4+7cnOY6pi6s39Aeuhk41XBG4L0UkQVUP5XH6x4ChwBjgd+BG3JKh+XkZaI77cL7Es+/vAl7zGvACrpDeY7hCcYtwCwVdCtwAvC4i36lqsud3uhxXJG4ScD+uWuqjuC9wI/O5Vh9gtapu8/Pc3cCnuFpI/YEnRaQebq2H24BYXOns1Z44AZ4BegB34pYhPQG3LCieOJvhSvcvAa7BFev7N24NmJNVVUWkCu7vpx7wEK5GVDPPltUPuBpSpiwo6QJattmWdcN/EcKs25Isxz6GKxQXk2XfyZ7jrvI8jvc87u95XBtX3G5UltdUwH2gaQGxPQVs8LN/Wo64rvNc89os+87z7Hs1y74awGHc2gv/397ZhGhVhXH892jEtDACh6yGsoSkRYkLixZt0iLChQUFLSWwUKmNKK1ibFeLdm0msEBEGgsb+tAGJoIM05iUUcIxYRTU7GMcP0iipvm3eM7N6/W+9700L71z5fnB4d5z3nPu+d/LC8895zn3POBG8RTwXuH6LybNCyu0DQO7SsrFtTurzsM3z5sqPLdB4EAufxR4paK/7cA4cHOu7H7gb2B1yr8MzADL2zzXtcA00NPt/1+k2aeYngq6wUX8LbiYPi3UewQYlnQpK5B0EDgJPNbi2g/h22gP5drM5PMdZCR3fiIdv8z1exEfrfSloqX4W/+gmd2UpdSmB3iwoq878EhxlTrSvU4Ao/nnlvT15fKHgc1mtsHMysIKPwHsBmZyOifwZ78i1VkJHJJ0uEI3Sfd8fMfYoOGE0Qi6wbQ84tk1CShOvdyJbyld5Gd8W+cysjn3XwrlxXwnuJA7/7OkLCvPfAm96fg5PgLJ0kQqv7uirx48xGg7HVmfVTrAo8x9DLwOjJvZj2b2Qu73Xjw88V+FtCSncyE+qmlHprunslbQCMKnEcxlfsLny4ssAkZbtDmXjrcD53PlZdf5v8n0vAQcKvl9oqQs3/a2TgmRdAH3d7xqZsuALcAOMxuT9EPqbzfu6ymSjXgmud5/UUam+3xlraARhNEI5jIHgPVmtkDSZQAzexj3Y+xr0eYIHnRoDe7HyJarrqnRX/FtvNOMA2fwQELvtqtc0va+zksCSWNmthkPzvQAvnhgBJ8uG5XU6qPIEeB5M1smaayii3uBSZU78YOGEUYjmMu8DawHvjCzN7m6euoI8FFZA0mTZjYAbDWzaTzm8brUth3HgEVmthZ3FP8m6eRsbyKnbcbMNgHbzexWYA9uqJYAzwDPSbrSovk3wLNmNi/5LWaFme3DRxJHcWf6OjyC28FUpT+df2Zm2/DRRR/wJPC+pK/wSHUbgeG0ZDkzbEslvZbrbgW+giq4AQifRjBnkfQr8Dg+ctgJvAN8jYfJbbXcFnyqZRs+X78TOIsboHYM4iul3sK/Lej/j9JbIukDfNSzHNiFf3exAfieq36RMoaAW/BYzp1gP76q6UP8vnuBp5WWBks6DjwKXAEGcAO3FfdPnEh1/sCd4Z8Ab6Q6W/DnDfy7tHoVLYx80Dwi3GsQNAQzGwJOS9rYbS11MbOncKN0l6Tfu60nmD1hNIKgISR/zgiwWNJUt/XUwcz2At9K6u+2lqAzxPRUEDQESd/h0z/3dFtLHdIX4/upNzUYNIQYaQRBEAS1iZFGEARBUJswGkEQBEFtwmgEQRAEtQmjEQRBENQmjEYQBEFQm38APXy1aV/LHl0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loss rate from exponential fit = 3.36829155446 (1/sec)\n"
     ]
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize = (6, 4))\n",
    "fileName_list = ['1288G', '980G']\n",
    "color_list = ['blue', 'red', 'black', 'darkgreen', 'orangered']\n",
    "\n",
    "working_dir = os.getcwd()\n",
    "callThis = join(working_dir, '1288G_molonly' + '.txt')\n",
    "out = np.loadtxt(callThis, delimiter='\\t', unpack=True)\n",
    "time = np.sort(out[0])\n",
    "num = out[9]\n",
    "num = num[np.argsort(out[0])]\n",
    "init_num = num[0]\n",
    "\n",
    "################\n",
    "## two-body fit\n",
    "popt2, pcov2 = curve_fit(twobody, time, num, p0=[20000, 1e-4])\n",
    "print(\"two body loss rate = {} (1/sec)\".format(popt2[-1]*1000.))\n",
    "twobody_coeff = popt2[-1]/popt2[0] ## this is beta constant divided by volume\n",
    "\n",
    "ax.scatter(time, num/init_num, c = 'k', label='1288G (NaLi only)')\n",
    "plt.plot(time, twobody(time, 1., popt2[-1]), c = 'k', linestyle='--')\n",
    "\n",
    "\n",
    "###################################\n",
    "def onetwobody(t, initnum, gamma1):\n",
    "    N0 = initnum\n",
    "    gamma2 = twobody_coeff \n",
    "    sol = odeint(fullModel, N0, t, args=(gamma1, gamma2))  \n",
    "    return sol[:, 0]\n",
    "###################################\n",
    "\n",
    "i = 0\n",
    "for fileName, color in zip(fileName_list, color_list):\n",
    "#     out = np.loadtxt(callThis, delimiter='\\t', unpack=True)\n",
    "    out = processData(join(working_dir, fileName + '.txt'))\n",
    "    \n",
    "    x = out[0]\n",
    "    num = out[1]\n",
    "\n",
    "    init_num = num[0]\n",
    "    num = num/init_num\n",
    "    \n",
    "    if flgerr == True:\n",
    "        err = out[2]\n",
    "        err = err/init_num\n",
    "    t = x\n",
    "    ax.scatter(t, num, c = color, label=fileName_list[i])\n",
    "    if flgerr == True:\n",
    "        ax.errorbar(t, num, err, c = color, ls='')\n",
    "    \n",
    "    #######\n",
    "    ## curve_fit to the full model with ode_int\n",
    "    popt, pcov = curve_fit(onetwobody, t, out[1], p0=[out[1][0], 1e-1])\n",
    "    \n",
    "    tt = np.linspace(0, 400, 100)\n",
    "    plt.plot(tt, onetwobody(tt, *popt)/popt[0], c = color)\n",
    "    print(\"for \" + str(fileName) + \", loss rate = {} (1/sec)\".format(popt[-1]*1000))\n",
    "    i += 1 \n",
    "    \n",
    "ax.set_ylabel('NaLi number',fontsize = 15)\n",
    "ax.set_xlabel('Hold time (msec)',fontsize = 15)\n",
    "ax.set_ylim([0., 1.2])\n",
    "ax.set_xlim([-5, 410])\n",
    "ax.legend()\n",
    "\n",
    "plt.xticks(fontsize = 12)\n",
    "plt.yticks(fontsize = 12)\n",
    "\n",
    "plt.savefig('lifetime comparison.pdf', dpi = 500)\n",
    "plt.savefig('lifetime comparison.png', dpi = 500)\n",
    "plt.show()\n",
    "\n",
    "########\n",
    "## simple one-body fit\n",
    "out = processData(join(working_dir, '1288G.txt'))\n",
    "popt, pcov = curve_fit(onebody, out[0], out[1], p0 = [out[1][0], 1e-2])\n",
    "print(\"loss rate from exponential fit = {} (1/sec)\".format(1000 * popt[-1]))\n",
    "\n",
    "# final = np.vstack((x_temp, num_temp, err_temp)).T\n",
    "# np.savetxt('aggregated.txt', final, delimiter ='\\t', fmt ='%f')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
